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ABSTRACT 

We present a new method to construct fully self-consistent equilibrium models of multi- 
component disc galaxies similar to the Milky Way. We define distribution functions for 
the stellar disc and dark halo that depend on phase space position only through action 
coordinates. We then use an iterative approach to find the corresponding gravitational 
potential. We study the adiabatic response of the initially spherical dark halo to the 
introduction of the baryonic component and find that the halo flattens in its inner 
regions with final minor-major axis ratios q = 0.75 - 0.95. The extent of the flattening 
depends on the velocity structure of the halo particles with radially biased models 
exhibiting a stronger response. In this latter case, which is according to cosmological 
simulations the most likely one, the new density structure resembles a “dark disc” 
superimposed on a spherical halo. We discuss the implications of these results for our 
recent estimate of the local dark matter density. 

The velocity distribution of the dark-matter particles near the Sun is very non- 
Gaussian. All three principal velocity dispersions are boosted as the halo contracts, 
and at low velocities a plateau develops in the distribution of v z . For models similar to 
a state-of-the-art Galaxy model we find velocity dispersions around 180 km s -1 for v z 
and the tangential velocity, v and 150 - 205kms _1 for the in-plane radial velocity, 
vr, depending on the anisotropy of the model. 
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1 INTRODUCTION 

It is now generally accepted that most of the rest mass in the 
Universe is contained in a type of “dark matter” that is inca¬ 
pable of electromagnetic interactions. Dark matter began to 
cluster gravitationally earlier than baryonic matter because 
it was decoupled from the cosmic radiation bath, which at 
redshifts > 1000 provided abundant pressure. Hence the for¬ 
mation of a galaxy starts with the accumulation of dark 
matter to form the dark halo. When and how baryons ac¬ 
cumulated in the gravitational potential wells of dark halos 
is ill understood. In the simplest picture, they simply fell in 
dissipatively and formed stars near each halo’s centre, but 
this picture is inconsistent with the small fraction of baryons 
that are now in stars rather than intergalactic gas. Several 
lines of evidence indicate that the release of nuclear energy 
by early stars strongly attenuated the infall of baryons, thus 
increasing the extent to which galaxies are dark-matter (dm) 
dominated. 

The least dm dominated galaxies are ones with lumi¬ 
nosities similar to that, L*, of our Galaxy: dwarf galaxies 
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are very strongly dm dominated and stars play such a sub¬ 
sidiary role in the most massive dm halos, those of galaxy 
groups and clusters, that we do not even recognise these 
halos as galaxies. Yet even in galaxies such as ours, it may 
be argued that the baryons have played a very subordinate 
role. Measurements of the optical depth to microlensing of 
stars located less than ~ 10 deg from the Galactic centre im¬ 
ply that a significant majority of the mass inside ~ 3kpc is 
concentrated in stars (Binney & Evans 2001; Bissantz et al. 
2004). But the majority of the gravitational force that holds 
the Sun in its orbit comes from dark matter (e.g. Piffl et al. 
2014, hereafter P14), and the dominance of dark matter in¬ 
creases sharply with increasing distance from the centre. 

Our understanding of how baryons accumulated at the 
centres of dark halos to form visible galaxies is limited. This 
is because the physics of this process is complex and is known 
to involve scales < 1 pc that are by far unresolved in feasi¬ 
ble simulations of the formation of galaxies in a cosmological 
context. Early simulations produced galaxies that contained 
too large a fraction of the baryons in the universe, and were 
much more dominated by their spheroidal components than 
is our Galaxy or our near neighbour M33. In recent years 
it has become clear that energy released during star forma- 
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tion is very efficient at preventing the accumulation of large 
masses of gas and stars in the centres of dark halos, by ex¬ 
pelling gas that has fallen in to a halo and the products 
of stellar evolution into intergalactic space, where the great 
majority of the baryons must reside. The high efficiency of 
such “feedback” from star formation is surprising and it can 
be reproduced in simulations only by modifying the encoded 
laws in essentially arbitrary ways. 

Since we cannot trust current simulations of the be¬ 
haviour of baryons during star formation, it is worthwhile 
to explore simple limiting cases. One such limiting case is 
that in which the baryons have accumulated in the Galaxy 
gradually rather than in lumps. This scenario is favoured by 
four observations. The first is that the bulge of our Galaxy 
bears all the hallmarks of having formed through a thin disc 
of stars that formed on nearly circular orbits in the Galactic 
plane developing a strong bar, which then buckled (Combes 
& Sanders 1981; Raha et al. 1991; Ness et al. 2013). The 
second relevant observation is that in the last ~ 10 Gyr the 
rate of star formation in the disc has declined by no more 
than a factor ~ 3 (e.g. Aumer & Binney 2009). A third rel¬ 
evant observation is that the star-formation rate in a disc 
appears to be proportional to the surface density of cold, 
molecular gas (Leroy et al. 2008). A fourth observation is 
that currently < 10% of the baryons in the disc are in gas. 

From the indication that all but a tiny fraction of the 
Galaxy’s stars formed from quiescently orbiting gas in the 
Galactic plane, we infer that our Galaxy has not experienced 
a significant merger since it started seriously forming stars. 
From the sustained rate of star formation over the Galaxy’s 
life together with the dependence of star formation rate on 
cold-gas density and the relatively small current stock of 
gas, we infer that the baryonic mass of our Galaxy has ac¬ 
cumulated over its entire life, presumably through sustained 
cooling of gas from the extended corona into the disc. There 
is even strong circumstantial evidence of this cooling process 
taking place now (Marasco et al. 2013). 

A dark halo will respond to the quiescent accumula¬ 
tion of baryons in a disc by distorting adiabatically: that 
is, the orbit of each dm particle will evolve in such a way 
that its action integrals J r , and J z will be constant, with 
the consequence that the density of dm particles in action 
space, /(J), is unaffected by the accumulation of baryons 
and the resulting evolution of the galaxy’s gravitational po¬ 
tential 3>(x) (e.g. Binney & Tremaine 2008, §4.6.1). The con¬ 
stancy of /(J) is enormously useful because the gravitational 
clustering of dark matter in isolation can be simulated with 
confidence, so we know what dark halos would now populate 
the Universe if there were no baryons. That is, from dm only 
simulations we can infer what /(J) would be in the absence 
of baryons, and the adiabatic principle assures us that /(J) 
will be unchanged by the quiescent infall of baryons in the 
observed quantities. Hence the hypothesis that the baryons 
accumulated quiescently so the dark matter responded adi¬ 
abatically to their arrival, allows us to predict the present 
structures of galaxies before we understand the complex pro¬ 
cesses that resulted in the accumulation of baryons. 

This general idea was appreciated at an early stage in 
the development of the cold dark matter (CDM) cosmolog¬ 
ical paradigm (Blumenthal et al. 1986), but until recently 
we have lacked the technology required to exploit it fully. 
The key to its exploitation is the ability to compute the 


action vector as a function J(x,v) of ordinary phase-space 
coordinates (x, v), for then it is possible to compute the dm 
density 

Pdm(x) = ^d 3 v/[J(x, v)] (1) 

at any location x. 

In an axisymmetric potential $ one action, = L z , is 
simply the component of angular momentum parallel to the 
potential’s symmetry axis, and in a spherical potential J z = 
L — 11/ 2 1, where L is the length of the angular-momentum 
vector L. Blumenthal et al. (1986) estimated the response 
of a dark halo to the accumulation of baryons by assuming 
that all dm particles are on circular orbits, so J r = 0. A 
more realistic estimate of halo response was obtained by 
Sellwood & McGaugh (2005), who lifted the restriction J r = 
0 to circular orbits but retained the assumption of spherical 
symmetry. The latter is awkward because the basic picture 
is of baryons accumulating in a disc, which will immediately 
break any pre-existing spherical symmetry of the dark halo. 

In this paper we lift the assumption of spherical sym¬ 
metry, so we are able to compute how the quiescent intro¬ 
duction of the disc would have deformed an initially spher¬ 
ical dark halo into a more centrally concentrated, oblate 
structure. We are, moreover, able to predict the full, three- 
dimensional distribution of the velocities of dm particles at 
any point in the Galaxy, in particular at the Sun. These 
distributions are triaxial and give non-Maxwellians speed 
distributions. 


2 METHODOLOGY 

Our strategy is to adopt dfs /(J) for the Galaxy’s stel¬ 
lar disc and dark halo and to solve for the gravitational 
potential 4?(x) that these self-consistently generate in the 
presence of given density distributions for the bulge and gas 
disc. We now describe our algorithm for the determination 
of $(x), which an extension of the methodology used by 
(Binney 2014, hereafter B14) to compute the observables of 
a flattened isochrone from an adopted df /(J). 

2.1 Iterative scheme 

B14 started from a guess $o at the self-consistent poten¬ 
tial. Using this guess and equation (1) he determined the 
density on a grid in the meridional ( R , z) plane, and then 
solved Poisson’s equation for the corresponding potential 
^ 1 / 2 (R,z). From this he made an improved guess 

$i{R, z) = (1 + a)$i/ 2 (R, z) - a$o {R, z) (2) 

at the self-consistent potential, where a <0.5 is parameter 
chosen to maximise the speed of convergence of the sequence 
of potentials to the self-consistent potential. 

The procedure just described involves a grid of points in 
the Rz plane at which the density is evaluated by integrating 
the df over v. Since the spatial scales of the disc and the 
dark halo are very discrepant, it is not cost-effective to use 
the same grid for both components: the grid would have 
to be fine enough to resolve the <0.1kpc structure of the 
disc and extensive enough to cover the dark halo, which has 
a significant amount of mass > 100 kpc from the Galactic 
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centre. Hence for each component 7 we use a different grid, 
and on this grid we store estimates of the density p 7 of that 
component and the contribution <f > 7 that it makes to the 
total gravitational potential 1 

fctot (R,Z)= <M R,Z )• (3) 

components -y 

Evaluation of the density of a single component involves 
evaluation of the actions J(x, v) and in this evaluation we 
use $ tot and the “Stackel Fudge” of Binney (2012a). Hav¬ 
ing evaluated the density p 7 of a component throughout the 
component’s grid, we solve Poisson’s equation for the cor¬ 
responding potential <f> 7 . The values of $ 7 at an arbitrary 
point (R, z) is obtained by interpolation on the values stored 
on the component’s grid. 

We start by determining the self-consistent density and 
potential implied by the adopted dark-halo df when the 
halo is isolated. Then our first trial potential for the whole 
Galaxy, , is the sum of this potential and the potentials 
of the bulge, gas disc and double-exponential stellar disc. We 
now compute the density of the dark halo in <f> t ot , and solve 
for the corresponding potential $dm and use this to update 
‘f’tot- Then we compute the density p* of the stellar disc in 
this estimate of and use the corresponding potential 4?* 
to update Then we compute the density of the dark 
halo in the updated ^tot, and so on until successive estimates 
of <3>tot differ negligibly. This occurs after no more than four 
computations of the density of each component. 

2.2 Reference model 

P14 modelled our Galaxy by combining a spheroidal mass 
distribution that has the classic NFW radial profile (Navarro 
et al. 1996) and a similar spheroidal bulge with a fixed gas 
disc and a dynamical stellar disc that has a df of the type 
introduced by Binney (2012b). The parameters of these com¬ 
ponents were adjusted until the composite model reproduced 
a variety of observational constraints, such as tangent veloc¬ 
ities uios (l) extracted from gas kinematics at various longi¬ 
tudes, the proper motion of Sgr A* and maser sources in 
the disc, and the kinematics of ~ 200 000 stars with spectra 
taken in the RAdial Velocity Experiment (RAVE Steinmetz 
et al. 2006; Kordopatis et al. 2013). Crucially, the stellar disc 
was also required to be consistent with the density profile 
above the Sun determined by Gilmore & Reid (1983) and 
Juric et al. (2008). 

P14 found that the most important unconstrained pa¬ 
rameter in their models was the axis ratio q < 1 of the 
dark halo. The smaller q is, the greater the fraction of the 
force on the Sun that comes from dark matter rather than 
stars. Since the potentials of the discs cause a dark halo that 
is spherical when isolated to flatten in its inner region, we 
have focused on the model in P14 that has a mildly flat¬ 
tened dark halo (minor-major axis ratio q = 0.8). Moreover, 
P14 remark that comparison of their results with those of 
Bienayme et al. (2014) favours a flattened model with axis 
ratio q ~ 0.8. We shall refer to this model as our ‘reference 
model’. It was described by P14 (cf. their figures 9, 10 and 

1 We actually store coefficients of Legendre-polynomial expan¬ 
sions and store also radial derivatives. 


11), but they did not specify its parameters. Appendix A 
gives the functional forms used for the density distributions 
of the model’s components. Table 4 contains the values of 
the parameters that appear in these forms. 

Many previous authors have, like P14, assumed that 
the dark halo will have the NFW radial profile. But this is 
the profile of dark halos in dm only simulations. Here we 
explore the extent to which the dark halo of the reference 
model would be modified if the discs determined by P14 
were adiabatically inserted into it. 

3 THE DFS 

3.1 The stellar disc 

The functional form of the stellar disc’s df, which is made up 
of contributions from the thick disc and each coeval cohort 
of thin-disc stars, was given by P14 and is reproduced in 
Appendix B. Table 5 gives the values of the parameters that 
appear in the disc’s df. 

We did make one minor change to the functional form 
of the DF: we imposed a lower limit, 1 kpc, on the circular ra¬ 
dius 7?c(J<)i) at which the epicycle frequencies are evaluated. 
Implicitly this introduces an upper limit on the frequen¬ 
cies k(J^), etc., that appear in the df. Crucially, after the 
limit has been imposed these functions remain well-defined 
functions of the actions, so the legitimacy of the df is not 
undermined. But on account of the limit, the frequencies em¬ 
ployed in the df differ from the epicycle frequencies at low 
J 4 ,. Imposing the limit has negligible impact on the struc¬ 
ture of the disc near the Sun; its material impact is at small 
radii, where the steep rise in the epicycle frequencies as the 
Galactic centre is approached would otherwise depress the 
disc’s density in an unphysical way. Table 5 gives the values 
of the parameters in the disc df. 

3.2 The stellar halo 

P14 used a df /( J) also for the stellar halo, but we omit this 
component because it has negligible mass and served only 
to improve fits to the stellar kinematics at small which 
we do not re-examine here. 

3.3 The dark halo 

Pontzen & Governato (2013) have fitted dfs to dark halos 
formed in dm only cosmological simulations. Unfortunately, 
the dfs they fitted are not usable here because they have 
the form f(E,L ): a df that describes one component of a 
composite model with a self-consistently generated poten¬ 
tial cannot have energy in its argument list for the following 
reason. When components are added, the central potential 
becomes deeper, so the energy of the star that rests at a com¬ 
ponent’s centre decreases. If E appears in the argument list, 
the phase space density around this orbit, J = 0, changes 
in an uncontrolled way. If / increases with |E| as it gen¬ 
erally does, the component’s density and mass increases, so 
\E\ increases, and without explicit mass renormalisation, the 
component’s mass is likely to run away to infinity during suc¬ 
cessive re-evaluations of the potential. When, by contrast, 
/ is taken to a function of just the actions, the phase-space 
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Figure 1. Lower panel: Density profiles computed from the 
dark halo DFs evaluated in the analytic NFW potential (coloured 
lines). The associated analytic NFW density profile is shown as 
a black dashed line. Upper panel: density profiles of the isolated 
halo models after the potential has been updated to the self- 
consistently generated one. 


density around each orbit never varies and the iterated po¬ 
tentials converge rapidly. 

In dm only simulations dark halos have NFW density 
profiles (Navarro et al. 1996). Posti et al. (2015) found a 
simple df /(J) that self-consistently generates a spherical 
model that has almost exactly an NFW profile. This model 
is essentially isotropic near its centre, and becomes mildly 
radially biased beyond its scale radius. For the dark halo 
we adopt dfs that are based on the form proposed by Posti 
et al. (2015) even though these dfs do not generate an NFW 
profile in the presence of the stellar disc. The normalisation 
of the dark-halo df (and thus the halo’s total mass) is chosen 
to match the circular speed, v C irc{Ro), in the reference model 
at the solar radius Ro- 

We now specify three different dark-halo dfs that cover 
the possibilities that the velocity distribution of dm particles 
is isotropic or radially or tangentially biased. We write 

r/fN _ N [1 + Jo/(dcore + h(J))]“ 

n ’ j 0 3 [i + fc(j)/j„F ’ 

where IV is a normalisation constant, a and f5 are con¬ 
stants that define the exponents of inner and outer power- 
law slopes of the density profile, Jo is a constant that con¬ 
trols the transition between the two. The constant J co re is 
a small number required to keep the central density finite. 
Finally, 

h(J) = Ij r + i^(|j„| + j 2 ) (5) 


Table 1 . DF parameters for the dark halo. 


Parameter 

Radially biased 

Isotropic 

Tangentially biased 

a 

5/3 

5/3 

5/3 

p 

2.9 

2.9 

2.9 

Jo 

47.3 

155 

218 

Jcore 

0.00025 

0.0001 

0.0002 

b 

8 

1 

0.125 


is an (approximately) homogeneous function 2 of J of order 
unity with 

6, + (6-i»)t»nh 2 (xTjnJTx) ' (6) 

B . s „ + (1-Mt»h-( X M±T_) (7) 

and bo = (1 + b)/2. In Appendix C we explain why realistic 
velocity distributions are obtained only when the values A 
and B are functions of the actions. 

The parameter b controls the model’s velocity 
anisotropy. For b = 1 we obtain a near-isotropic model, 
while values below (above) unity yield tangentially (radially) 
biased models. The quantities k and fare epicycle and 
azimuthal frequencies in the self-consistent spherical model 
evaluated at the radius R c of a circular orbit with angular 
momentum 


Jtot = Jr + |Jip| + Jz- (8) 

We take the argument of R c to be Jtot to make it an approx¬ 
imate function of energy, so R c does not become small, and 
the epicycle frequencies large, for stars on eccentric and/or 
highly inclined orbits with small \J V \. 

Following Posti et al. (2015), we set the exponents a 
and /3 to 5/3 and 2.9, respectively, to obtain a good ap¬ 
proximation to the NFW density profile. Table 1 lists all 
the adopted values of the parameters that define the halo’s 
df. They generate a halo that (a) closely approximates the 
NFW profile when the halo is isolated, and (b) in the pres¬ 
ence of the discs and bulge has a flattening interior to the 
Sun similar to that of the reference model’s dark halo. 

Once we have obtained a self-consistent model of an iso¬ 
lated halo, we freeze the dependence of R c on its argument, 
so while we relax <I>tot onto the potential that is jointly gen¬ 
erated by all the Galaxy’s components, the df stays exactly 
the same function of J. It is essential to freeze the function 
/(J) during the introduction of the discs and the bulge if one 
seeks to learn how the halo is distorted by the gravitational 
fields of its companions. 


3.4 The bulge and gas disc 

We could adopt a df /(J) for the final Galactic component, 
the bulge, but we do not do so for two reasons. First, the 
bulge contains a rotating bar and at present we cannot rep¬ 
resent such an object with a df /(J) (but see Sanders & 
Binney 2015a, for a non-rotating triaxial example). Second, 


2 Posti et al. (2015) argued that h(J) should be a homogeneous 
function of J, i.e. they required that h(aj) = a h ( J ). 
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Figure 2. Velocity anisotropy of the DM particles as a function 
of radius in the equatorial plane. 



7 = d In p / d In r 



Figure 4. Upper panel: solid thin lines show the spherically sym¬ 
metric density profiles of halos in isolation, while thick lines show 
the density profiles in the rry-plane after the halos have been con¬ 
tracted by introducing the baryons. The dotted lines show the 
profiles obtained by contracting the original models with the clas¬ 
sic adiabatic prescription. Lower panel: final density divided by 
initial density when the contraction is done properly (full curves) 
and using the classical prescription (dotted curves). 


Figure 3. Density slope - velocity anisotropy relation for our 
radially biased (purple) and isotropic models (orange) as well as 
predictions from cosmological dark matter-only simulations. 

we are primarily interested in the structure of the dark halo 
in the vicinity of the Sun, which should be well recovered 
using an axisymmetrised form of the bulge. Hence we repre¬ 
sent the bulge by the fixed density distribution specified in 
Appendix A. We also represent the Galaxy’s gas disc by a 
fixed density distribution. 

4 RESULTS 
4.1 Isolated halos 

We consider three halo dfs, which differ only in the value 
of the parameter 6 : we set 6 = 0.125 to generate a tangen¬ 
tially biased model, b = 1 for a nearly isotropic model, and 
6 = 8 for a radially biased model. In all plots we distin¬ 
guish these models through the colour scheme 6=1 orange 
isotropic, 6 = 8 green radially biased, and 6 = 0.125 pur¬ 
ple tangentially biased. The lower panel of Fig. 1 shows the 
density profiles generated by the three dfs in the analytic 
NFW potential (which has no core). The df parameters were 
optimized to provide a good match between these density 
profiles and that of the NFW halo fitted by P14. The up¬ 
per panel of Fig. 1 shows the radial density profiles of these 


models after the potential has relaxed to self-consistency. 
(Fitting a self-consistent DF-potential pair to the analytic 
profile directly is conceptually straightforward, but compu¬ 
tationally expensive, and is beyond the scope of this paper.) 

The thin curves in Fig. 2 show the anisotropy parameter 

/3a = 1 “ hjf) (9) 

for the initial isolated halo as a function of radius in the xy 
plane. To provide some context, Fig. 3 shows /3 a as a function 
of the slope of the logarithmic density profile for two of these 
models. The grey band shows the range of such dependency 
found by Ludlow et al. (2011) from cosmological simulations, 
while the black line shows the linear dependence that Hansen 
& Moore (2006) found described cosmologically simulated 
halos well. We see that our isotropic and radially biased 
models bracket the empirical results well, except possibly 
at the smallest values of |-y |, which occur at the centre of a 
halo. 


4.2 Halo distortion 

The introduction of the baryonic components causes the ini¬ 
tially spherical dark halo to contract. The end point of this 
contraction should coincide perfectly with what we would 
have found had we slowly grown the baryonic components 
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Figure 5. Minor-major axis ratio of the halo component as a 
function of elliptical radius a = K 2 + (z/q)' 2 . 


in a live N-body simulation of the halo. The almost coin¬ 
cident heavy curves in the upper panel of Fig. 4 show the 
final, relaxed DM density profiles in the disc (xy) plane. Even 
though the central parts of the initial profiles differ, with the 
radially-biased model being most cuspy (light full curves in 
the upper panel), the final profiles are very similar. The full 
curves in the lower panel of Fig. 4 show the ratio of initial 
to final density at each radius. We see that in all three cases 
(radial bias, isotropic, tangential bias), the central density 
increases by more than a factor 8. 

This figure shows (i) that at r < 3 kpc the baryonic com¬ 
ponent distorts the dark halo substantially, and (ii) that the 
variation of the initial density profiles with the anisotropy 
parameter b is just what is required to ensure that the three 
profiles coincide after the baryons have been added: as Sell- 
wood & McGaugh (2005) showed, the tangentially biased 
halo is more strongly compressed than the radially-biased 
halo, so before the barons are added its density needs to be 
lower than that of the radially-biased halo. 

The dotted curves in the upper panel of Fig. 4 show 
the profiles one obtains when one starts with the three iso¬ 
lated halos and applies to these objects the classic adiabatic 
contraction formalism of Blumenthal et al. (1986). This for¬ 
malism is based on the fiction that all dark matter parti¬ 
cles are on circular orbits, i.e. an extreme case of tangential 
bias. From the lower panel of Fig. 4 we see that the ana¬ 
lytic prediction is in excellent agreement with our results for 
r > 1 kpc, but under-estimates the contraction for smaller 
radii. 

We can for the first time study how the flattened baryon 
distribution breaks the spherical symmetry of the halo, and 
pulls it towards the plane hitherto analytic prescriptions 
for adiabatic distortion have been restricted to the spherical 
case. Fig. 5 shows the axis ratio of the final dark halo as a 
function of elliptical radius. We see that the radially biased 
halo is most strongly flattened by the baryons - its axis 
ratio does not rise above 0.8 until near Ro . The axis ratio 
of the initially nearly isotropic halo rises more gradually 
from similar values at ~ 0.2 kpc to ~ 0.9 near the Sun. We 
emphasise that these axis ratios are arising naturally and 
are not related to our choice of a reference model with a 
similarly flattened halo. In fact, we found very similar axis 


ratios in test runs using the P14 model with a spherical halo 
and then inspired by these results chose to use the flattened 
model as our reference. 

Fig. 6 shows the fractional change in the dark halo’s 
density (pfl na i — pini) / pini that is induced by adiabatic inser¬ 
tion of the baryons. We see that for any velocity anisotropy 
the effect is strongly concentrated to the inner Galaxy, 
r < 5 kpc. In the case of radial velocity anisotropy, the dark 
halo’s density is most enhanced in the region |o|<lkpc, 
and the enhancement could be (mis-)interpreted as the for¬ 
mation of a dark disc around the baryonic disc. 

4.3 Rotation curve 

Our P14 reference model was fitted to data constraining the 
Galaxy’s rotation curve using an NFW halo profile. As we 
have seen above, the adiabatic contraction of the dark halo 
significantly alters the halo’s central density profile. To coun¬ 
teract these changes a far as possible, we re-scaled the halo’s 
mass such that we still obtained the right circular speed in 
the solar cylinder. However, because contraction steepens 
the halo’s inner density profile, the halo’s contribution to 
rises inwards faster than in the P14 model, and when the 
contracted halo is used in conjunction with the P14 disc and 
bulge the model no longer matches the constraints on our 
Galaxy’s inner rotation curve. 

Comparison of the black and full green curves in Fig. 7 
shows the extent of the mismatch. In the left panel, the 
model circular-speed curve remains flat down to much lower 
radii than the curve derived from the observations. The right 
panel shows the extent of the conflict with observations of 
the terminal velocities of HI gas (Malhotra 1995). This is a 
fundamental problem that ties back to the use of an NFW 
profile in our reference model. This profile fits baryon-free 
dark halos and will not fit dark halos to which baryons have 
been adiabatically added. If the addition of baryons is non- 
adiabatic, the NFW profile might fit real dark halos, but we 
have no reason to suppose that the addition of baryons was 
non-adiabatic to the precise extent required for the NFW 
profile to remain valid. 


4.4 DM velocity distributions 

We now examine the velocity distribution of the dark mat¬ 
ter particles within the Galactic plane of our models. Fig. 2 
shows the radial dependence of the classical anisotropy pa¬ 
rameter (equation 9) before and after adiabatic contraction. 
The most significant changes are in the innermost regions 
where the baryonic discs are dominant. In the radially bi¬ 
ased case (green curves), the anisotropy increases slightly at 
R < 15 kpc and scarcely changes further out. In the other 
two cases the anisotropy increases at all radii, but no quali¬ 
tative changes are observable outside the disc region. 

Fig. 8 shows the distributions of the three velocity com¬ 
ponents (after marginalising over the other two components) 
of dm particles at the location x = (iio,0) of the Sun. 3 The 


3 We adopt Rq = 8.3 kpc. The distributions at the actual solar 
position 14 to 25 pc above the Galactic plane are indistinguishable 
from these. 
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Figure 6. Relative density difference of the dark halo component before and after the adiabatic contraction, (pfl na i — pini)/pini - due 
to the baryonic components. Each panel corresponds to a model with different velocity anisotropy as indicated above the panels. The 
contour lines mark locations where the density has increased similarly by a given fraction in per cent as indicated on the lines. 




Figure 7. Left: Circular speed curve of the target mass distribution (black line) and our final self-consistent model with a radially biased 
velocity structure in the dark halo. The other two models yield almost indistinguishable curves. The individual contributions from the 
dark halo and the baryonic component are also shown. Right: Terminal velocity curve as a function of Galactic longitude, /, predicted 
by the same models as in the left panel and for comparison the HI terminal velocity measurements from Malhotra (1995) to which the 
P14 models were fitted. 


Table 2. Local dark matter velocity dispersions 


bias 

o- VR 

kms -1 

Cvtp 

kms -1 

&VZ 

kms -1 

radial 

174 

150 

154 

isotropic 

154 

154 

157 

tangential 

141 

156 

158 


non-Gaussian nature of these velocity distributions is evi¬ 
dent. In all cases adiabatic contraction of the halo broadens 
the distributions for the same reason gas that is adiabati- 
cally compressed heats up. As the indirection broadens, its 
top becomes remarkably flat. The dispersions of the distri¬ 
butions are listed in Table 2. 

The sensitivity of particle detectors typically depends 
on the speed of the particles to be detected (e.h. Green 
2012). The full curves in Fig. 9 show for each final model 
the distribution of the speeds of dm particles with respect to 
the Sun. The tangentially-biased model has the most slow- 
moving and least fast-moving particles, while the radially 



Heliocentric speed v [km s x ] 

Figure 9. Full curves: the distributions speeds with respect to 
the Sun of DM particles at the location of the Sun in our three final 
models. Dashed curves show the Maxwellian distribution with the 
same dispersions. 
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vr [fans *] v v [kms - 1 ] v z [kms x ] 


Figure 8. Dark matter velocity distribution in the solar annulus for our three halo models. 


Table 3. Parameters of the best-fitting Tsallis speed distributions 
shown in Fig. 10 

bias q vq J k m s ~ 1 ] 

radial 0.834 341.8 

isotropic 0.837 332.7 

tangential 0.825 331.8 


biased model has ~ 40% more particles with v > 600 kms -1 . 
Hence dm particles will be more readily detected in the 
radially-biased case, which is fortunately the case to which 
N-body simulations point (e.g. Hansen & Moore 2006; Lud¬ 
low et al. 2011; Pontzen & Governato 2013). 

The dashed curves in Fig. 9 shows the Maxwellian dis¬ 
tributions with the same velocity dispersions as the real dis¬ 
tributions. These Maxwellian equivalents are quite similar to 
one another and differ materially from the actual distribu¬ 
tions. In any stellar system the actual velocity distribution 
has to fall below a Maxwellian at large v because it has to 
vanish above the escape speed, and this effect is evident in 
Fig. 9. 

A model distribution of speeds that falls to zero at a 
finite speed is that proposed by Tsallis (1988): 


n(v) oc v 2 


i - (i -q) : 


9 /( 1 - 9 ) 


( 10 ) 



Figure 10. The full curves are the same as in Fig. 9. The dashed 
curves show fits to the speed distribution provided by the Tsallis 
function (10). Table 3 gives the parameters of the models. 


invisible under the true distribution. The fit to the tangen¬ 
tially biased model is very good, but the Tsallis distribution 
completely fails to reproduce the distribution of the radially- 
biased model. 


4.5 The stellar disc 


which tends to a Maxwellian as q —¥ 1. Ling et al. (2010) 
showed that the Tsallis distribution fitted the speeds of par¬ 
ticles in a spherical shell of radius r ~ 8 kpc in their cosmo¬ 
logical hydrodynamical simulations of baryon infall. The full 
curves in Fig. 10 are identical to those in Fig. 9, while the 
dashed curves show the best-fitting Tsallis distributions. Ta¬ 
ble 3 gives the parameters of these fits. The fit to the speed 
distribution of the isotropic (orange) model is excellent - 
in fact, the curve for the Tsallis distribution is frequently 


Figs. 11 and 12 show the stellar density distribution in the 
case of the radially biased halo by showing the disc’s density 
at fixed z as a function of R and at Rq as a function of z, 
respectively. The black curves show the result of integrating 
the disc df in the potential of the reference model, while the 
green curves are obtained by integrating the same df in the 
model with the contracted halo.. The differences between 
these profiles are significant only in the outer disc (Fig. 11) 
where the final halo is rounder than that of the reference 
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Figure 11. Density profile of the stellar disc in slices of constant 
z before and after the dark halo’s response to the presence of the 
disc. The line pairs correspond to heights above the Galactic plane 
from 0 kpc (up-most almost identical lines) to 10 kpc (lowest 
lines) in steps of 2 kpc. 



Figure 12. Vertical density profile of the stellar disc in the solar 
cylinder before and after the response of the halo to the presence 
of the disc. The star count data from Juric et al. (2008) have 
been re-normalised (shifted vertically) to allow a comparison to 
the model. 

halo. Fig. 12 shows that in the solar cylinder the two density 
profiles agree extremely well, so the final model still provides 
an excellent match to the star count data of Juric et al. 
(2008) to which the reference model was fitted. 


5 RELATION TO OTHER WORK 

The growth of baryonic discs in dark halos has been ex¬ 
tensively studied with particle simulations (e.g. Pillcpich 
et al. 2014, and references therein). In the great majority 
of studies, the manner in which baryons have infiltrated the 
dark halo has been controlled by “sub-grid physics”, that is, 
prescriptions governing gas cooling, star formation, and the 
consequent heating and ejection of gas. These prescriptions 
are only loosely based in physical principles and are freely 
adjusted to bring certain statistical properties of the bary¬ 
onic structures that form into agreement with observations. 
In these simulations dark halos do not necessarily evolve 
adiabatically - in fact in most simulations it is likely that 


scattering of dm particles by bars and massive gas clouds 
will have significantly lowered the peak phase-space density 
of the dark halo. 

Debattista et al. (2008) adiabatically introduced rigid 
discs into triaxial dark halos that had previously been 
formed by colliding spherical dark halos. After their adi¬ 
abatic introduction, the discs were adiabatically removed 
by evaporation, and the final distributions of dark-matter 
particles were compared with the initial distribution. The 
principal conclusion of this study was that the initial and 
final distributions were very similar, which indicates that 
the dark halos responded adiabatically. In the absence of a 
disc, the halos were mostly strongly prolate, and their outer 
regions remained strongly prolate at all times, but the discs 
distorted the inner regions of their halos to near axisymmet- 
ric oblate forms. Thus Debattista et al. (2008) underlines the 
value of being able to compute the response of a dark halo to 
the adiabatic insertion of axisymmetric baryonic structures. 

Another study in which for a period the dark halo must 
have evolved adiabatically is that of DeBuhr et al. (2012), for 
between redshift 2 = 1.3 and 2=1 they adiabatically grew 
a massive disc. At 2 = 1 the disc became a fully fledged N- 
body disc capable of developing a bar and strict adiabaticity 
ended. The paper focuses on the axis ratios of the halo and 
the evolution of the disc. It provides no quantitative infor¬ 
mation regarding the evolution of the phase-space density 
of dm particles. 

Pillepich et al. (2014) compared the final distribution 
of dark matter in two simulations of the formation of an 
object that is slightly less massive than our Galaxy, in one 
of which baryons were included while the other contained 
only dark matter. The radial DM profiles for the models 
with and without baryons in their Figure 1 are very simi¬ 
lar to our plots of the DM density before and after baryon 
addition in the upper panel of Fig. 4. They found that the 
inclusion of baryons increased dm density at the Sun by ~ 30 
per cent, in excellent agreement with our findings. Two pro¬ 
cesses contribute to enhancement of the dm density in the 
plane of the simulation with baryons. One is the pinching 
of the dark halo by the disc’s gravitational field, and the 
other is the capture of satellites onto highly inclined orbits 
through dynamical friction on the baryonic disc (Quinn & 
Goodman 1986) and their subsequent tidal shredding into a 
rotating dark disc (Abadi et al. 2003; Read et al. 2008, 2009; 
Ruchti et al. 2014; Read 2014). Our models include the first 
mechanism but not the second. The second process forms 
a dark disc that rotates at a significant speed in the same 
sense as the baryonic disc, while pinching of the dark halo 
does not set the halo rotating. Pillepich et al. (2014) showed 
that for the rather quiescent merger history expected for our 
Galaxy, the effect of pinching predominates, accounting for 
~ 2/3 of the overall dark disc. 

Ling et al. (2010) and Pillcpich et al. (2014) used hy- 
drodynamical simulations of baryon insertion to examine 
the velocity distribution of dm particles at the Sun’s loca¬ 
tion. The simulation of Pillepich et al. (2014) involved 13 
million DM and 13 million baryonic particles in the high- 
resolution region, while the simulations of Ling et al. (2010) 
employ slightly fewer particles. Pillepich et al. (2014) in¬ 
ferred the local velocity distribution of DM particles from 
the ~ 80 000 DM particles in the region |i? — Ro\ < 2kpc, 
1 2 1 < 2 kpc, while Ling et al. (2010) studied the 2 662 par- 
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tides in the smaller volume |i? — i?o| < lkpc, \z\ < 1 kpc. 
Despite the enormous computational resources deployed to 
run these simulations, in neither study is Poisson noise unim¬ 
portant, and Pillepich et al. (2014) considered it necessary 
to smooth their velocity distribution to a 50kms _1 resolu¬ 
tion. Our models require a tiny fraction of the computational 
resource and provide velocity distributions that relate pre¬ 
cisely to the Sun’s location and yet are completely free of 
Poisson noise. Neither Ling et al. (2010) nor Pillepich et al. 
(2014) were able to study the dependence of halo kinematics 
on the extent of radial bias in the initial dark halo. 

The speed distributions at the Sun reported by Ling 
et al. (2010) and Pillepich et al. (2014) are not unlike the 
one we obtain in the initially isotropic case. A significant 
difference between their results and ours is that they find 
significant net rotation of the local DM particles whereas 
we have none. For example Ling et al. (2010) find ( v</,} ~ 
35kms _1 . Streaming in the dark halo is a consequence of 
non-adiabatic interaction between the baryons and the dark 
halo. In addition to the shredding of dark halos discussed 
above, a bar will set a dark halo spinning by losing angular 
momentum to it (Tremaine & Weinberg 1984). 


6 DISCUSSION & CONCLUSIONS 

We have shown how to construct fully self-consistent multi- 
component galaxy models with axisymmetric action-based 
distribution functions /(J). As a worked example we have 
taken a Galaxy model presented by P14 and made its dark 
halo a fully dynamical object. In particular, we determined 
the gravitational potential that is self-consistently generated 
by the stellar disc and the dark halo in the presence of pre¬ 
defined contributions from a gas disc and an axisymmetric 
bulge. 

To make the dark halo a dynamical object, we must se¬ 
lect its velocity anisotropy. To this end we have introduced 
a new family of dfs. We have used the new dfs to explore 
three options: radially biased, isotropic and tangentially bi¬ 
ased dark halos. The first two models bracket the predictions 
of cosmological DMonly simulations. 

The most important change in the final model is the 
adiabatic compression of the dark halo by the flattened po¬ 
tential of the discs; we can study this process at a tiny frac¬ 
tion of the computational cost of a cosmological simulation. 
The contraction has two components. The first component 
is a spherical shrinkage. At radii r > 1 kpc, i.e., roughly out¬ 
side one disc scale height, this is quite well approximated by 
the simple formalism of Blumenthal et al. (1986). At smaller 
radii we find stronger contraction. This is the reverse of what 
is usually found with N-body simulations, in which generally 
less contraction is observed than Blumenthal et al. predict 
(Abadi et al. 2010; Gnedin et al. 2010). 

The second component of compression is a pinching to¬ 
wards the plane at radii r < 10 kpc, which flattens the halo. 
The extent to which the halo flattens on introduction of the 
discs depends of its velocity anisotropy in the sense that ra¬ 
dial bias maximises the flattening. Our initially spherical, 
radially biased model deforms to an axis ratio q < 0.8 out 
to ~ 5 kpc. The tangentially biased halo, by contrast, has 
q > 0.9 at r > 0.3 kpc. This has important consequences, 
since cosmological simulations predict radially biased con¬ 


figurations. Indeed, the left panel of Fig. 6 shows that in the 
case of a radially biased dark halo, the difference between 
the actual, compressed halo and the original spherical halo 
looks very much like a “dark disc”. A “disc” of this type 
would be non-rotating, in contrast to a dark disc formed by 
the capture of satellite halos onto low-inclination orbits fol¬ 
lowed by tidal shredding (Quinn & Goodman 1986; Abadi 
et al. 2003; Read et al. 2008; Ruchti et al. 2014; Read 2014). 

Our models yield the detailed velocity distribution of 
DM particles at the Sun, upon which depends the detectabil¬ 
ity of dark matter in direct detection experiments. The dis¬ 
tributions of speed with respect to the Sun are distinctly 
non-Gaussian, being more and less sharply peaked in the 
radially and tangentially biased models, respectively. Given 
the tendency of particle detectors to have a threshold speed 
for detection, the radially biased model offers significantly 
better chances of detecting dark matter. 

Recently, P14 modelled the Galaxy by combining dy¬ 
namical models of the stellar disc with star counts and the 
kinematics of RAVE stars to obtain estimates the local DM 
density. They found that the major model uncertainty was 
the flattening of the dark halo. The best-fitting local dm den¬ 
sity increased with the assumed halo axis ratio q as oc q~ 0 ' 89 
(cf. their Fig. 9). Here we find that even if the dark halo were 
originally spherical, the dark halo must be flattened inside 
the solar radius Ro, with q = 0.7 to 0.9. Even though q 
varies with radius, this variation is modest within Ro, the 
relevant region from the perspective of P14. 

Using a very different approach to that of P14, Bien- 
ayme et al. (2014) determined the local dm density from a 
subset of the data used by P14. The local dm densities of 
these two studies agree for halo axis ratio q in the range 0.79 
to 0.94. It is interesting that this is just the range of axis 
ratios the present study leads us to expect if the dark halo 
were spherical before the disc was added. 

The dfs of our dark halos were chosen so they gener¬ 
ate isolated dark halos that closely follow the NFW profile. 
When cohabiting with the P14 disc, the halos becomes more 
centrally concentrated, with the consequence that the final 
circular-speed curves do not match the data as well as in the 
reference model, where the dark halo has an NFW profile in 
the presence of the discs and bulge. The models do, however, 
have appropriate values v c (Ro) for the circular speed at the 
Sun. 

There are two distinct ways in which we could modify 
the present model so as to restore a satisfactory fit to the 
circular-speed curve while keeping the dark halo alive. One 
way is to compensate for the increased central concentra¬ 
tion of the halo by increasing the disc scale lengths from the 
values adopted by P14. The other is to modify the df of 
the dark halo in such a way that it approximates the NFW 
profile when its density is evaluated in the presence of the 
baryons rather than when isolated. The case for modifica¬ 
tion of the df is that the introduction of the baryons was 
not completely adiabatic, with the consequence that /(J) 
decreased near the origin of action space and increased away 
from the origin to leave the total halo mass (2-7t) 3 f d 3 J /(J) 
constant. The modification of the df could be used to in¬ 
troduce a degree of rotational streaming such as will arise 
from friction against the bar and entrapment of sub-halos 
on highly inclined orbits. Cosmological simulations could be 
used to guide the choice of halo df. 
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While there is a strong case for modifying the halo df 
from the very simple form adopted here, it is not clear that 
the modified df should yield an NFW profile in the presence 
of the baryons. Therefore, the next step should be to see 
whether a plausible disc df exists which when combined 
with the present halo df will produce a satisfactory model 
of our Galaxy. 

We hope soon to report new constraints on the struc¬ 
ture of our Galaxy’s dark halo obtained by fitting to ob¬ 
servational data model Galaxies in which the stellar disc, 
stellar halo and dark halo are all represented by dfs of the 
form /(J). The time is also ripe to fit such models to data 
for external galaxies, which a new generation of integral-field 
spectrographs are making extremely rich. We anticipate that 
these models will constrain the dm distributions strongly be¬ 
cause the stars will be described by either distinct dfs for 
each observationally distinguishable range of metallicities, 
or by an “extended distribution function” (Sanders & Bin- 
ney 2015b). 
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APPENDIX A: THE MASS MODEL 

Our mass models have five components: a gas disc, a thin 
disc, a thick disc, a flattened bulge and a dark halo. Since 
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the stellar halo has negligible mass, we do not include it in 
the mass model. For the density laws of the disc components 
we have 


p{R, z) = exp 


R \z\ -Rhole \ 
Rd z d R ) 


( 11 ) 


A non-zero parameter Rhoie creates a central cavity in the 
disc. This is used to model the gas disc while for the other 
two discs .Rhoie is set to zero. The values of the parameters 
are given in Table 4. 

The density distributions of the dark halo and the bulge 
components are 

P(jR ’ Z) = rrU(l + m )f>- 7 ex PI - ( mr o/ r cut) 2 ], (12) 


where 


m(R, z) = \J(R/ro) 2 + ( z/qr 0 ) 2 . (13) 


The parameters of the reference model are given in Table 4. 


APPENDIX B: THE DISC DF 

The disc DF is built up out of “quasi-isothermal” compo¬ 
nents. The DF of such a component is 

/(Jr, j„ L z ) = f ar (Jr, L Z )U Z (J z , L z ), (14) 

where f tTr and f az are defined to be 

nr 2 

fa r (J r ,L z ) = — 5 —[1 + tanh(L z /Lo)]e _reJr/CTr (15) 

7r O r K, 

and 

f aM (J„L,) = ^-ze- vJ ‘rt. (16) 

Here Lt(L z ), k{L z ) and v{L z ) are, respectively, the circular, 
radial and vertical epicycle frequencies of the circular orbit 
with angular momentum L z , while 

S(L 2 ) = E 0 e- flc/iid , (17) 

where R C (L Z ) is the radius of the circular orbit, determines 
the surface density of the disc: to a moderate approxima¬ 
tion this surface density can be obtained by using for L z in 
equation (17) the angular momentum L Z (R) of the circular 
orbit with radius R. The functions a r {L z ) and a z {L z ) con¬ 
trol the radial and vertical velocity dispersions in the disc 
and are approximately equal to them at R c . Given that the 
scale heights of galactic discs do not vary strongly with ra¬ 
dius (van der Kruit & Searle 1981), these quantities must 
increase inwards. We adopt the dependence on L z 

a r (L z ) = aroe (R °- Rc)/R "’ r 

a z (L z )=a z0 e (R °- R ' )/R (18) 

so the radial scale-lengths on which the velocity dispersions 
decline are R a ,i- Our expectation is that I? CT y ~ 2 Rd- 

In equation (15) the factor containing tanh serves to 
eliminate retrograde stars; the value of La controls the ra¬ 
dius within which significant numbers of retrograde stars 
are found, and should be no larger than the circular angular 
momentum at the half-light radius of the bulge. Provided 
this condition is satisfied, the results for the extended solar 
neighbourhood presented here are essentially independent of 

To- 


Table 5. Parameters for stellar disc DF. 




Thin disc 

Thick disc 

cr r 

kms -1 

33.9 

50.5 


kms -1 

25.0 

48.9 

Rd 

kpc 

2.58 

2.58 

R(r,r 

kpc 

9 

12.9 

Rtr,z 

kpc 

9 

4.1 

Tl 

Gyr 

0.01 

- 

Tm 

Gyr 

10 

- 

TO 

Gyr 

8 

- 

P 

- 

0.33 

- 

-^thick 

- 

- 

0.455 


We take the DF of the thick disc to be a single pseudo- 
isothermal. The thin disc is treated as a superposition of 
the cohorts of stars that have age r for ages that vary from 
zero up to the age r max — 10 Gyr of the thin disc. We take 
the DF of each such cohort to be a pseudo-isothermal with 
velocity-dispersion parameters oy and oy that depend on age 
as well as on L z . Specifically, from Aumer & Binney (2009) 
we adopt 


oy -{L z , t) — oyo 


( t + n 

V T m + Tl 


g(^0 — Re) jR(T,r 


oy (L z ,r) — (7 z 0 


( t + n 

V Tm + Tl 


e (i?o— Rc) l Ra ,z 


(19) 


Here cr 2 o is the approximate vertical velocity dispersion of 
local stars at age r m ~ 10 Gyr, n sets velocity dispersion 
at birth, and (3 ~ 0.33 is an index that determines how 
the velocity dispersions grow with age. We further assume 
that the star-formation rate in the thin disc has decreased 
exponentially with time, with characteristic time scale to, so 
the thin-disc DF is 


/thn(Jr, Jz,L z ) 


J 0 Tm dr e^ f ° f„ r (J r , L z )f am (J z , L z ) 
to(e Tm R° — 1 ) 


( 20 ) 

where ay and a z depend on L z and r through equation (19). 
We set the normalising constant Eo that appears in equation 
(17) to be the same for both discs and use for the complete 

DF 


/disc {Jr, Jz, L z ) — fthn{Jr, Jz, Lz') T Fthk/thk ( Jr , J z , Lz'), 

( 21 ) 

where F t hk is a parameter that controls the fraction (1 + 
^thk) — 1 °f stars that belong to the thick disc . 4 The values 
of the parameters are given in Table 5. 


APPENDIX C: CUSPS AND DIMPLES IN 14 
DISTRIBUTIONS 

To implement velocity anisotropy in our halo DF we did not 
opt for the simplest possible solution. Here we explain the 
reasons for this decision. In the DF the two factors A and 


4 Note, that Fthk is the ratio of the total masses of the thick and 
the thin discs, while the parameter /thk used for the mass model 
is the ratio of the local mass densities of the two discs. Hence the 
two parameters are intimately related but not the same. 
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Figure 13. Local v<p distributions using a simpler form of the 
halo DF with a constant anisotropy parameter. 


B (Eqs. 6, 7) control the anisotropy of the system by en¬ 
hancing or suppressing the population of orbits with certain 
properties. These factors are constant except for very radial 
orbits. This variation makes h( J) not a strictly homogeneous 
function any more, a property that was demanded by Posti 
et al. (2015) when they proposed this class of dfs. 

If we remove the J dependency by substituting the tanli 
function with unity we obtain A = b and B = 1 which 
is sufficient to create an anisotropic model. However, when 
we examine the tangential velocity distribution of such a 
model (Fig. 13) we find that in the radially biased model, 
the distribution in v v has a cusp at v$ = 0, while in the 
tangentially biased model it has a dimple at v$ = 0; only the 
isotropic model has a v$ distribution of the expected domed 
shape. The cusp and dimple become more pronounced after 
adiabatic contraction. 

The cusp and the dimple arise because the energy, upon 
which the ergodic DF of the isotropic model depends, is a 
quadratic function of v#, while J# = Rv^ is a linear func¬ 
tion of Vt/,. So when the DF of the isotropic model is written 
as a function of the actions and then Taylor expanded in v$ 
at fixed x, the linear dependence of / upon v$ that arises 
from the dependence of / on has to be cancelled by linear 
dependencies of J r and J z on v</,. If these linear dependen¬ 
cies do not cancel, df /dv$ is non-zero at v$ = 0, and P V(j> 
acquires a cusp or a dimple there. Below we show that this 
delicate cancellation is disturbed by shifting b from unity. 

The derivative of the v v distribution P Vip plotted in 
Fig. 8 is 


dP„, 

dv<p 




II 


AVrAV; 


df_ 

dvu, 


( 22 ) 


and 


df_ 

dV<n 


df dJ r df dJ v df dJ z 
dJ r dv v dJ v dv v dJ z dv v 


(23) 


Now 


and 


df_ 

dJr 


df 1 


dh(J) b 


df_ = df_ = df n v 
dJ„ dJ z dh{ J) K 


(24) 


(25) 


We require derivatives with respect to at a fixed location 
x = (R, 0) in the plane, so 


dJ v 

dvu, 


= R, 


(26) 


where we have used J v = Rv v . 

We don’t have analytic expressions for dJ r /dv v and 
dJ z /dv v , but the following argument shows that dJ r /dv v < 
0 for 0 < Vip Vdrc(R). Consider a star with low v v , and 
therefore a small guiding radius, that is just able to reach 
(R, 0), i.e. has its apocentre there. It is on an eccentric orbit 
with substantial J r . Increasing v v moves the guiding centre 
outwards and thus diminishes the amplitude of the star’s 
radial oscillations and lowers its value of J r . 

To see that dJ z /dv v | x < 0 for 0 < v v <C v c \rc{R), 
consider the case of a spherical potential. Then J z = L— \ J ( j > \ 

with L = R Jv^ + uf in the Galactic plane, so 


dJ z 
dv$ 


^L- R = ^±- R = R f'I± 

oVcf, L \L 



(27) 


which is clearly negative for small v^. In a strongly flattened 
potential, the vertical motion largely decouples from the pla¬ 
nar motion, so dJ z /dv$ becomes numerically small, but it 
remains negative. 

Inserting all these insights into equation (23) we arrive 
at 


df df fidJr n y dj z \ 

dv v dh{ J) \b dv v k k dv v ) 

for the integrand in equation (22). Since df /dh( J) < 0, the 
sign of AP Vv /Avtp is opposite to the sign of the bracket in 
equation (28). In the isotropic case (b = 1) we must have 
AP Vv ,/Av v ~ 0 at v v = 0, so 


d Jr + Sl<p dJ z ^ Slip R 
dv v k dv v K 


(29) 


Hence for 6 ^ 1 this balance is disturbed with the conse¬ 
quence that P v</> acquires a non-zero gradient at v# = 0, the 
gradient being negative when b > 1 and positive when b < 1 
as is evident in Fig. 13. 

While cusps and dimples are not strictly unphysical, 
we think it unlikely that real dm distributions would exhibit 
them. A requirement to avoid these features amounts to a 
restriction on permissible forms of the function h(J) as, e.g. 
the one we choose for this work. 


© 2014 RAS, MNRAS 000, 1 -13 















